function[IRF_VAR_boot_org,IRF_LP_boot_org,IRF_VAR_boot,IRF_LP_boot,IRF_VAR_boot_full,IRF_LP_boot_full]=VAR_bootstrap_replacement(Y_var,hor,p,identification,n_boot,responseV,normalizeV)



[~,rhs,lhs]=lag_variable(Y_var,p);
[T,N]=size(lhs);
[M]=size(rhs,2);

phi=(rhs'*rhs)\(rhs'*lhs);
e_org=lhs-rhs*phi;

% Center the residuals (as Lütkepohl recommends)
e_mean = mean(e_org);                   % sample mean
e = e_org - e_mean;            % centered residuals
% e=e_org;

for i_boot = 1:n_boot  %1000 bootstraps 
    
    lhs_boot = NaN(T,N);
    rhs_boot = NaN(T,M);

    %Set initial conditions
    init      = unidrnd(T);
    rhs_boot(1,:) = rhs(init,:);
    
    for ii = 1:T
        % Get a random draw of the empirical residuals
        lhs_boot(ii,:) = phi'*rhs_boot(ii,:)' + e(unidrnd(T),:)';
        if ii < T
            rhs_boot(ii+1,:) = [1 lhs_boot(ii,:) rhs_boot(ii,2:end-N)];
        end
    end
  

    %Estimate IRFs for VAR and LP 
   [IRF_VAR_boot_full(:,:,i_boot),~,~,~]=VAR_model(lhs_boot,hor,p,identification,0);

   [IRF_LP_boot_full(:,:,i_boot)]=LP_Jorda_model(lhs_boot,hor,p,identification,0);

end

%Org and Scale bootstrap IRFs
for i_boot=1:n_boot
    IRF_VAR_boot(i_boot,:) = IRF_VAR_boot_full(responseV,:,i_boot) / IRF_VAR_boot_full(normalizeV,1,i_boot);  
    IRF_LP_boot(i_boot,:) = IRF_LP_boot_full(responseV,:,i_boot) / IRF_LP_boot_full(normalizeV,1,i_boot); 

    IRF_VAR_boot_org(i_boot,:) = IRF_VAR_boot_full(responseV,:,i_boot);
    IRF_LP_boot_org(i_boot,:) = IRF_LP_boot_full(responseV,:,i_boot);

end



end

